This lecture gives an overview of some optimization tools in Julia.
1.1 Flowchart
Statisticians do optimizations in daily life: maximum likelihood estimation, machine learning, …
Category of optimization problems:
Problems with analytical solutions: least squares, principle component analysis, canonical correlation analysis, …
Problems subject to Disciplined Convex Programming (DCP): linear programming (LP), quadratic programming (QP), second-order cone programming (SOCP), semidefinite programming (SDP), and geometric programming (GP).
Nonlinear programming (NLP): Newton type algorithms, Fisher scoring algorithm, EM algorithm, MM algorithms.
Large scale optimization: ADMM, SGD, …
Flowchart
1.2 Modeling tools and solvers
Getting familiar with good optimization softwares broadens the scope and scale of problems we are able to solve in statistics. Following table lists some of the best optimization softwares.
LP
MILP
SOCP
MISOCP
SDP
GP
NLP
MINLP
R
Matlab
Julia
Python
Cost
modeling tools
cvx
x
x
x
x
x
x
x
x
x
A
Convex.jl
x
x
x
x
x
x
O
Optimization.jl
x
x
x
x
x
x
x
O
JuMP.jl
x
x
x
x
x
x
x
O
MathOptInterface.jl
x
x
x
x
x
x
x
O
convex solvers
Mosek
x
x
x
x
x
x
x
x
x
x
x
A
Gurobi
x
x
x
x
x
x
x
x
A
CPLEX
x
x
x
x
x
x
x
x
A
SCS
x
x
x
x
x
x
O
COSMO.jl
x
x
x
x
O
Hypatia.jl (more cones)
x
x
x
x
O
NLP solvers
NLopt
x
x
x
x
x
x
O
Ipopt
x
x
x
x
x
x
O
KNITRO
x
x
x
x
x
x
x
x
$
O: open source, A: free academic license, $: commercial
Modeling tools such as cvx (for Matlab) and Convex.jl (Julia analog of cvx) implement the disciplined convex programming (DCP) paradigm proposed by Grant and Boyd (2008) http://stanford.edu/~boyd/papers/disc_cvx_prog.html. DCP prescribes a set of simple rules from which users can construct convex optimization problems easily.
Solvers (Mosek, Gurobi, Cplex, SCS, COSMO, Hypatia, …) are concrete software implementation of optimization algorithms. My favorite ones are: Mosek/Gurobi/SCS for DCP and Ipopt/NLopt for nonlinear programming. Mosek and Gurobi are commercial software but free for academic use. SCS/Ipopt/NLopt are open source.
Modeling tools usually have the capability to use a variety of solvers. But modeling tools are solver agnostic so users do not have to worry about specific solver interface.
If you want to install the commercial solvers Gurobi or Mosek, instructions are below:
Gurobi: 1. Download Gurobi at link. 2. Register an account and request free academic license at link. 3. Run grbgetkey XXXXXXXXX command on terminal as suggested. It’ll retrieve a license file and put it under a specified folder, e.g., /Users/huazhou/Documents/Gurobi/gurobi.lic. 4. Set up the environmental variables. On my machine, I put following two lines in the ~/.julia/config/startup.jl file: ENV["GUROBI_HOME"] = "/Library/gurobi1101/macos_universal2" and ENV["GRB_LICENSE_FILE"] = "/Users/huazhou/Documents/Gurobi/gurobi.lic".
Mosek: 1. Request free academic license at link. The license file will be sent to your edu email within minutes. Check Spam folder if necessary. 2. Put the license file at the default location ~/mosek/.
Install Julia packages Convex.jl, SCS.jl, Gurobi.jl, Mosek.jl, MathProgBase.jl, NLopt.jl, Ipopt.jl, which are open source.
1.3 DCP Using Convex.jl
Standard convex problem classes like LP (linear programming), QP (quadratic programming), SOCP (second-order cone programming), SDP (semidefinite programming), and GP (geometric programming), are becoming a technology.
DCP Hierarchy
1.3.1 Example: microbiome regression analysis
We illustrate optimization tools in Julia using microbiome analysis as an example.
16S microbiome sequencing techonology generates sequence counts of various organisms (OTUs, operational taxonomic units) in samples.
Microbiome Data
For statistical analysis, counts are normalized into proportions for each sample, resulting in a covariate matrix \(\mathbf{X}\) with all rows summing to 1. For identifiability, we need to add a sum-to-zero constraint to the regression cofficients. In other words, we need to solve a constrained least squares problem \[
\text{minimize} \frac{1}{2} \|\mathbf{y} - \mathbf{X} \beta\|_2^2
\] subject to the constraint \(\sum_{j=1}^p \beta_j = 0\). For simplicity we ignore intercept and non-OTU covariates in this presentation.
Let’s first generate an artifical data set.
usingRandom, LinearAlgebra, SparseArraysRandom.seed!(257) # seedn, p =100, 50X =rand(n, p)# scale each row of X sum to 1lmul!(Diagonal(1./vec(sum(X, dims=2))), X)# true β is a sparse vector with about 10% non-zero entriesβ =sprandn(p, 0.1) y = X * β +randn(n);
1.3.2 Sum-to-zero regression
The sum-to-zero contrained least squares is a standard quadratic programming (QP) problem so should be solved easily by any QP solver.
Suppose we want to know which organisms (OTU) are associated with the response. We can answer this question using a sum-to-zero contrained lasso \[
\text{minimize} \frac 12 \|\mathbf{y} - \mathbf{X} \beta\|_2^2 + \lambda \|\beta\|_1
\] subject to the constraint \(\sum_{j=1}^p \beta_j = 0\). Varying \(\lambda\) from small to large values will generate a solution path.
usingConvex# # Use Mosek solver# using Mosek# solver = Mosek.Optimizer# # Use Gurobi solver# using Gurobi# solver = Gurobi.Optimizer# # Use SCS solver# using SCS# solver = SCS.Optimizer# Use Hypatia solverusingHypatiasolver = Hypatia.Optimizer# # Use Clarabel solver# using Clarabel# solver = Clarabel.Optimizer# solve at a grid of λλgrid =0:0.01:0.35# holder for solution pathβ̂path =zeros(length(λgrid), size(X, 2)) # each row is β̂ at a λ# optimization variableβ̂classo =Variable(size(X, 2))# obtain solution path using warm start@timefor i in1:length(λgrid) λ = λgrid[i]# define optimization problem# objective problem =minimize(0.5sumsquares(y - X * β̂classo) + λ *sum(abs, β̂classo))# constraint problem.constraints +=sum(β̂classo) ==0# constraintsolve!(problem, solver; silent =true) β̂path[i, :] = β̂classo.valueend
┌ Warning: Concatenating collections of constraints together with `+` or `+=` to produce a new list of constraints is deprecated. Instead, use `vcat` to concatenate collections of constraints.
└ @ Convex /Users/huazhou/.julia/packages/Convex/y7lu0/src/deprecations.jl:129
1.819383 seconds (4.05 M allocations: 381.876 MiB, 2.31% gc time, 70.21% compilation time: 14% of which was recompilation)
usingCairoMakief =Figure()Axis( f[1, 1], title ="Sum-to-Zero Lasso", xlabel = L"\lambda", ylabel = L"\beta_j")series!(λgrid, β̂path', color =:glasbey_category10_n256)f
1.3.4 Sum-to-zero group lasso
Suppose we want to do variable selection not at the OTU level, but at the Phylum level. OTUs are clustered into various Phyla. We can answer this question using a sum-to-zero contrained group lasso \[
\text{minimize} \frac 12 \|\mathbf{y} - \mathbf{X} \beta\|_2^2 + \lambda \sum_j \|\mathbf{\beta}_j\|_2
\] subject to the constraint \(\sum_{j=1}^p \beta_j = 0\), where \(\mathbf{\beta}_j\) are regression coefficients corresponding to the \(j\)-th phylum. This is a second-order cone programming (SOCP) problem readily modeled by Convex.jl.
Let’s assume each 10 contiguous OTUs belong to one Phylum.
# # Use Mosek solver# using Mosek# solver = Mosek.Optimizer# # Use Gurobi solver# using Gurobi# solver = Gurobi.Optimizer# # Use SCS solver# using SCS# solver = SCS.Optimizer# Use Hypatia solverusingHypatiasolver = Hypatia.Optimizer# # Use Clarabel solver# using Clarabel# solver = Clarabel.Optimizer# solve at a grid of λλgrid =0.0:0.005:0.5β̂pathgrp =zeros(length(λgrid), size(X, 2)) # each row is β̂ at a λβ̂classo =Variable(size(X, 2))@timefor i in1:length(λgrid) λ = λgrid[i]# loss obj =0.5sumsquares(y - X * β̂classo)# group lasso penalty termfor j in1:(size(X, 2) ÷10) βj = β̂classo[(10(j -1) +1) :10j] obj = obj + λ *norm(βj)end problem =minimize(obj)# constraint problem.constraints = [sum(β̂classo) ==0] # constraintsolve!(problem, solver; silent =true) β̂pathgrp[i, :] = β̂classo.valueend
1.972572 seconds (5.89 M allocations: 588.133 MiB, 3.00% gc time, 62.61% compilation time: 51% of which was recompilation)
We see it took Mosek <2 second to solve this seemingly hard optimization problem at 101 different \(\lambda\) values.
usingCairoMakief =Figure()Axis( f[1, 1], title ="Sum-to-Zero Group Lasso", xlabel = L"\lambda", ylabel = L"\beta_j")series!(λgrid, β̂pathgrp', color =:glasbey_category10_n256)f
1.3.5 Example: matrix completion
Load the \(128 \times 128\) Lena picture with missing pixels.
We fill out the missin pixels uisng a matrix completion technique developed by Candes and Tao \[
\text{minimize } \|\mathbf{X}\|_*
\]\[
\text{subject to } x_{ij} = y_{ij} \text{ for all observed entries } (i, j).
\] Here \(\|\mathbf{M}\|_* = \sum_i \sigma_i(\mathbf{M})\) is the nuclear norm. In words we seek the matrix with minimal nuclear norm that agrees with the observed entries. This is a semidefinite programming (SDP) problem readily modeled by Convex.jl.
This example takes longer because of high dimensionality. COSMO.jl seems to be the fastest solver for this problem. Other solvers take excessively long time.
# Use COSMO solver (fast)usingCOSMOsolver = COSMO.Optimizer# # Use Hypatia solver (slow)# using Hypatia# solver = Hypatia.Optimizer()# # Use Clarabel solver (slow)# using Clarabel# solver = Clarabel.Optimizer()# Linear indices of obs. entriesobsidx =findall(Y[:] .≠0.0)# Create optimization variablesX =Variable(size(Y))# Set up optmization problemproblem =minimize(nuclearnorm(X))problem.constraints += X[obsidx] == Y[obsidx]# Solve the problem by calling solve@timesolve!(problem, solver) # fast
------------------------------------------------------------------
COSMO v0.8.9 - A Quadratic Objective Conic Solver
Michael Garstka
University of Oxford, 2017 - 2022
------------------------------------------------------------------
Problem: x ∈ R^{32897},
constraints: A ∈ R^{41025x32897} (41281 nnz),
matrix size to factor: 73922x73922,
Floating-point precision: Float64
Sets: DensePsdConeTriangle of dim: 32896 (256x256)
ZeroSet of dim: 8128
Nonnegatives of dim: 1
Settings: ϵ_abs = 1.0e-05, ϵ_rel = 1.0e-05,
ϵ_prim_inf = 1.0e-04, ϵ_dual_inf = 1.0e-04,
ρ = 0.1, σ = 1e-06, α = 1.6,
max_iter = 5000,
scaling iter = 10 (on),
check termination every 25 iter,
check infeasibility every 40 iter,
KKT system solver: QDLDL
Acc: Anderson Type2{QRDecomp},
Memory size = 15, RestartedMemory,
Safeguarded: true, tol: 2.0
Setup Time: 27.64ms
Iter: Objective: Primal Res: Dual Res: Rho:
1 -1.4585e+03 1.5985e+01 5.9854e-01 1.0000e-01
25 1.4525e+02 5.1105e-02 1.1356e-03 1.0000e-01
50 1.4758e+02 1.1725e-02 1.4834e-03 6.8658e-01
75 1.4797e+02 5.5160e-04 4.7489e-05 6.8658e-01
100 1.4797e+02 1.7304e-05 1.3870e-06 6.8658e-01
------------------------------------------------------------------
>>> Results
Status: Solved
Iterations: 100
Optimal objective: 148
Runtime: 1.147s (1147.0ms)
2.172648 seconds (3.98 M allocations: 717.390 MiB, 4.64% gc time, 62.44% compilation time: 15% of which was recompilation)
Problem statistics
problem is DCP : true
number of variables : 1 (16_384 scalar elements)
number of constraints : 1 (8_128 scalar elements)
number of coefficients : 8_128
number of atoms : 3
Solution summary
termination status : OPTIMAL
primal status : FEASIBLE_POINT
dual status : FEASIBLE_POINT
objective value : 147.9711
Expression graph
minimize
└─ nuclearnorm (convex; positive)
└─ 128×128 real variable (id: 104…399)
subject to
└─ == constraint (affine)
└─ + (affine; real)
├─ index (affine; real)
│ └─ …
└─ 8128×1 Matrix{Float64}
We use MLE of Gamma distribution to illustrate some rudiments of nonlinear programming (NLP) in Julia.
Let \(x_1,\ldots,x_m\) be a random sample from the gamma density \[
f(x) = \Gamma(\alpha)^{-1} \beta^{\alpha} x^{\alpha-1} e^{-\beta x}
\] on \((0,\infty)\). The loglikelihood function is \[
L(\alpha, \beta) = m [- \ln \Gamma(\alpha) + \alpha \ln \beta + (\alpha - 1)\overline{\ln x} - \beta \bar x],
\] where \(\overline{x} = \frac{1}{m} \sum_{i=1}^m x_i\) and \(\overline{\ln x} = \frac{1}{m} \sum_{i=1}^m \ln x_i\).
1.4.1 Define NLP optimization problem using Optimization.jl